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The mechanism of diffusion in supercooled liquids is in- 
vestigated from the potential energy landscape point of view, 
with emphasis on the crossover from high- to low-T dynamics 
over the range Ta > T > T c . Molecular dynamics simulations 
with a time dependent mapping to the associated local min- 
inum or inherent structure (IS) are performed on unit-density 
Lennard- Jones (LJ). New dynamical quantities introduced in- 
clude r2i a (t), the mean-square displacement (MSD) within a 
basin of attraction of an IS, R2(t), the MSD of the IS itself, 
and g(t), the distribution of IS waiting times. The config- 
uration space is treated as a composite of the contributions 
of cooperative local regions, and a method is given to obtain 
the physically meaningful gioc{t) and mean waiting time ri oc 
from g(t). A new understanding of the crossover is obtained in 
terms of r2; 3 (i) and n oc . At intermediate T, r2j s (t) posesses 
an interval of linear i-dependence allowing calculation of an 
intrabasin diffusion constant Di S . Near T c , where intrabasin 
diffusion is well established for t < n oc , diffusion is intrabasin 
dominated with D — Di S ; D may be calculated within a basin. 
Below T c , Ti oc exceeds the time, t p i, needed for the system 
to explore the basin, indicating the action of barriers at the 
border; r; oc = t v i is a criterion for a transition to activated 
hopping. Intrabasin diffusion provides a means of confine- 
ment not involving barriers and plays a key role in dynamics 
above T c . The distinction between motion among the IS (IS 
dynamics) below T c and saddle, or border dynamics above T c , 
where the system is always close to one of the saddle-barriers 
connecting the basins and IS boundaries are closely spaced 
and easily crossed, is discussed. A border index is introduced 
based upon the relation of R2(t) to the conventional MSD, 
and shown to vanish at T ~ T c . It is proposed that intra- 
basin diffusion is a manifestation of saddle dynamics. 



I. INTRODUCTION; DIFFUSION AND THE 
POTENTIAL ENERGY LANDSCAPE 

Knowledge of the physical mechanism of diffusion is 
essential to an understanding of the dynamical complex- 
ities of supercooled liquids. A mechanism should include 
the collective motions composing the elementary diffusive 
step, any correlations among the steps, and a statement 
about whether barrier crossing is important. If so, distri- 
butions of barrier heights, 'reaction coordinates' for bar- 
rier crossing, and details of the relevant saddle-barriers 
are significant. Valuable information is encoded in the 
T-dependence of the self diffusion coefficient D and of 
the viscosity rj (the Stokes-Einstein law, D ~ T/rj, holds 
until very close to the glass transition at T g ). For most 
normal liquids D exhibits Arrhenius T-dependence. In 
jjj strong liquids this continues all the way down to T g , 



while super- Arrhenius, with the apparent activation en- 
ergy Ea(T) increasing with decreasing T, sets in before 
T g in Jjj fragile liquids. Some exceptions classified as 
'weak' by Kivelson et. al H, including simple atomic 
models such as Lennard- Jones (LJ), show a normal or 
upper-supercooled range of power-law D(T). Clearly the 
mechanism of diffusion is T-depcndent, and a straightfor- 
ward conclusion is that activated potential-energy barrier 
crossing, or hopping, plays an important role. In con- 
trast to chemical reactions where the Arrhenius law was 
first introduced, however, the barriers and reaction co- 
ordinates for diffusion have proved elusive. To further 
complicate matters, free energy or 'entropic' barriers, as- 
sociated with the difficulty of finding a low-energy path- 
way, may also be important. It is not straightforward to 
infer a mechanism from D(T), and thus there has been 
considerable speculation. 

Adam and Gibbs suggested || that liquids may be 
decomposed into roughly independent cooperatively re- 
arranging regions (CRR) containing z* particles, with 
super- Arrhenius arising from the growth of the CRR, and 
thus the activation free energy, with decreasing T . An in- 
fluential paper by Goldstein Q formulated dynamics in 
terms of the SA^-dimensional potential energy [/-surface, 
or landscape. He proposed the existence of distinct high 
and low-T dynamical regimes, separated by a crossover 
temperature. At low T the system stays for times long 
compared to a vibrational period in the basins associated 
with the local minima of the landscape, and diffusion 
is governed by infrequent hopping of the saddle-barriers 
connecting adjoining basins. Perhaps the rearrangements 
of the CRR correspond to reaction coordinates on the 
landscape. At high-T a different, freer motion is pre- 
sumed to occur. The onset of the low-T mechanism is 
sometimes described as a 'transition to hopping', but if 
the high-T diffusion constant is Arrhenius, implying ac- 
tivation, the change must be more subtle. The idea of a 
crossover is central to current thinking about supercooled 
liquids. We will pursue the landscape point of view in this 
paper, so the mechanism of diffusion is naturally discused 
in terms of the topology. 

Goldstein's picture has been extended to involve two 
characteristic temperatures, meaning that the mecha- 
nism changes gradually over a crossover range Ta > T > 
T c . Interesting features of supercooled dynamics such 
as super-Arrhenius and stretched exponential decay of 
time correlations first appear at Ta Q and D(T) extrap- 
olates to zero from above, via power-law fit, at the mode- 
coupling H temperature T c < Ta- Note that T c , defined 
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by the extrapolation, seems a useful quantity regardless 
of the validity of mode-coupling theory. An alternative 
formulation, the 'frustration limited domain' theory of 
Kivelson and co-workers, introduces an upper crossover 
temperature T* at which collective dynamics first be- 
comes important. Still lower f^jj^ lie T g , the Kauzmann 
temperature Tk and the VTF To, but this paper is con- 
cerned with temperatures from Ta to somewhat under 
T c . It is difficult to obtain well-equilibrated computer 
simulation data below T c , and impossible to approach 

Stillinger and Weber showed |^|| how to obtain an ex- 
plicit description of motion over the landscape in a molec- 
ular dynamics simulation, by mapping a 3A-dimensional 
liquid configuration r(i) to the associated minimum R(t), 
or inherent structure (IS), with frequent steepest-descent 
minimizations (faster conjugate-gradient minimizations 
are equivalent [jlOj] at supercooled T). The configuration 
space is thus partitioned into the basins of attraction of 
the IS and the rate of inherent structure transitions (1ST) 
may be measured. The IS-mapping is a powerful tech- 
nique for the investigation of supercooled dynamics but it 
has not yet been exploited to anything like its full poten- 
tial. Sastry et al. demonstrated JllJ] in the LJ mixture 
that the averaged IS energy, < Ui S (T) >, undergoes a 
sharp drop (the || 'slope') from an upper high-T to a 
lower low-T plateau. They identify the beginning of the 
drop with Ta, and T c lies near the apparent bottom. The 
bottom energy is cooling-rate dependent and the true 
value, requiring cooling slow enough to maintain equili- 
bration, remains unknown. It is appealing, within the 
landscape point of view, to relate the change in dynami- 
cal mechanism to occupation by the system of lower-lying 
IS. 

More direct investigations of the influence of the land- 
scape on dynamics are possible. The minima of the 
squared potential gradient W = \WU\ 2 include all the ex- 
terna, or critical points, of U, including minima of order 
K = and saddles of order K > 0. Minimizing W defines 
a mapping to a configuration which need not be a min- 
imum of U. With the critical point mapping, Cavagna, 
Angelani et al. [|l2|,[l3| showed that minima dominate be- 
low T c , saddles above. The IS-mapping assigns the sys- 
tem to a minimum even if it is in the upper reaches of the 
basin near the saddle-barriers; the critical point mapping 
does so only if it is actually close to the IS. They suggest 
||12||l3[ that the crossover is from motion among the min- 
ima (IS dynamics) at T < T c to motion among the sad- 
dles (saddle dynamics) at T > T c . Either mapping may 
be used at any T. Reference to saddle dynamics or IS 
dynamics below indicates that it gives the simpler, more 
physical description. The saddles ruled |l2| regime is also 
one of border dynamics. Reaction coordinates leading to 
several different IS intersect at a multidimensional or j^] 
'monkey' saddle, so IS boundaries are closely spaced in 
the vicinity of a saddle. The system is always near a bor- 



der, crossings are facile and frequent, and 1ST are not 
hops but merely unphysical 'bookkeeping' jl0| events. A 
much better understanding of the T > T c region is being 
achieved with the new ideas about saddles and border 
dynamics. 

Recently |l0| we gave the first quantitative expres- 
sion for D in terms of the 1ST rate < u>i s > and the 
mean-square separation of successive IS. We argued that 
the long-time slope of the mean-square displacement 
per degree of freedom (MSD) of the IS-configuration, 
R2(t) =< (AR(i)) 2 > /6iV, is equal to D, as is so for the 
ordinary MSD, r2(i). The IS-displacement is the sum 
of the 1ST- vectors (separations of successive minima), 
AR(t) = ^Za=i SH a after n(t) transitions in time t; for 
any quantity x, the displacement Ax(t) = x(t) — x(0). A 
key quantity in this approach is the IST-vector cor- 
relation, C(/3) =< (<5R a • 5H. a +p)u)i S >. In unit-density 
LJ, N = 32 atoms, where T c = 0.52 (LJ units), an IST- 
Markov approximation of random walking among the IS, 
predicts D accurately for T <T C and substantially over- 
estimates D for T > T c . The Markov approximation 
is equivalent to the statement that the IST-vectors are 
uncorrelated, C(0) = 0,/3 > 0. Our result supports the 
idea of long sojourns in the basins below T c , since the sys- 
tem then has time to lose correlation or memory between 
1ST, leading to a random walk. Because of the uniquely 
weak T-dependence of D at T > T c in LJ, 'transition 
to hopping' may actually be an accurate description of 
the crossover in this case. The overestimate of D by the 
Markov approximation at T > T c is associated [jLOl with 
a long ranged anticorrelation of successive transition vec- 
tors, negative C(/3) at large /?, necessary for frequent 1ST 
to result in relatively little displacement. Anticorrelated 
IST-vectors are a signature of the 'bookkeeping' 1ST oc- 
curing in border dynamics. 

A problem arises in trying to learn about landscape 
dynamics with the IS-mapping. In the absence of long- 
range correlations the configuration space of the liquid 
should be regarded as a composite of the spaces of N crr 
independent local regions of z* ~ 0(1) particles. A 
simple model |lj] uses 3Af sinusoidal potentials. Local 
regions and Adam-Gibbs H CRR need not be equiv- 
alent, but we will treat them as such. With indepen- 
dent 1ST occurring in each region, and with a transition 
for the whole system recorded when any region changes, 
< LUi S >^ O(N) and the distribution of waiting times 
between 1ST, g(t), is correspondingly ||, and unphysi- 
cally, skewed towards short time ~ 0(1/N). Computa- 
tional constraints then indicate that one should simulate 
the smallest realistic system; with sufficiently large N an 
1ST will be observed on every time step and no mean- 
ingful calculation will be possible. More fundamentally, 
the desired properties are those of the CRR, the build- 
ing blocks of the macroscopic liquid, but straightforward 
simulation yields those of the homogenized composite. 
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With increasing N the former become more difficult to 
disentangle from the latter. This is also true for static 
quantities. For example, whatever the distribution of Ui S 
in a local region, it will be Gaussian for the entire sys- 
tem at large N . Even for the best choice of N some new 
theory is required to interpret the simulation data. 

In this article we further examine the T-dcpendent 
mechanism of diffusion in unit-density LJ, N = 32, via 
three new landscape dynamical entities, and with an ini- 
tial attempt to derive local relaxation from simulated 
composite dynamics. The first is r2i s (t) =< (Ar(t))f s > 
/6N, the MSD with no contribution from 1ST. Starting 
in a thermal ensemble of basins, r2^ s (t) is calculated from 
only those trajectories which remain in the initial basin 
at time t. The second is R2(t) and the third is the distri- 
bution gio C {t) of waiting times in an IS for a local region 
or CRR. We present a simple method to obtain gi oc (t) 
from g(t) and N crr , estimated from the averaged partic- 
ipation ratio, < Pr >, of the 1ST- vectors. The num- 
ber of degrees of freedom (not particles) participating 
in an 1ST is roughly < Pr >; thus z* s»< Pr > /3, 
N crr w 3 A/ < Pr>. 

Comparison of r2j S (i), R2(t), r2(t) and gi oc (t) casts 
considerable light upon the relation of diffusion to the 
landscape. For t >> t p i, where t p i is a characteristic 
plateau time, r2j s (i) attains a constant value, due to 
|L5| the finite volume of the basin. This is true even 
at high T where the system moves freely among the IS, 
since only trajectories which remain in the original IS 
contribute (thus averaging is difficult). At T = 1.10 the 
plateau is reached quickly but a small interval of lin- 
ear i-dependence is visible at intermediate times, from 
which an intrabasin diffusion coefficient, £)j s , may be 
calculated. Below Ta ~ 1-0 the diffusive interval in 
r2i S (t) expands; Ta is estimated as the top of the slope in 
< Ui S {T) >. More cooling reveals that the T ~ T c regime 
is intrabasin dominated with D — Di S - the true diffusion 
constant may be calculated in a single basin! Thus, as for 
high T, diffusion cannot consist of a sequence of activated 
1ST hops, but the mechanism is quite different from the 
high-T mechanism. Intrabasin diffusion arises from the 
roughness of the landscape at the upper elevations of a 
basin and provides a mechanism for confinement differ- 
ent from inter-basin barriers but nonetheless capable of 
producing |l(| a Markov chain of 1ST under the right 
conditions. We suggest that the elementary steps in this 
process are transitions among the domains of the sad- 
dles connected to the IS, ie, intrabasin diffusion is saddle 
dynamics. Intrabasin dominated diffusion is also IS dy- 
namics, since the motion among the IS is then a simple 
random walk; at intermediate T both descriptions are 
appropriate. 

The system does not sense the finite size of the basin 
until the mean local waiting time, r; oc , reaches t p i- We 
thus have a new criterion for the onset of the low-T IS- 
dynamics mechanism and, indeed, it is found that n oc 



crosses t p i slightly below T c . The analysis would make 
no sense with an A-dependent mean time taken from 
g(t). More information is obtained by comparison of 
R2(t) with r2(t). Let r(i) = R(t) + u(t), where Q 
u(i) is the return vector, the separation of the system 
from its IS. If AR and Au are uncorrelated, r2(t) must 
lie above R2(t). This is so at low T but as T increases 
R2(t) rises above r2(t). The relation of R2(t) to r2(t), in- 
dicating the correlation of AR and Au, yields a border in- 
dex B(T) which complements the information contained 
in C{0) regarding IS vs saddle dynamics. In sum, with 
some novel applications of the IS- mapping we achieve an 
extremely detailed description of the T-dependent mech- 
anism of diffusion. 

II. SIMULATION METHODS 

Molecular dynamics simulations are performed on 
unit-density supercooled L J |Io| , [l6| , |l7| , A = 32, with the 
methods described in ref |10|. As discussed there, and 
above, a small system should be used for 1ST dynamics. 
Increasing N beyond the number required to 'solvate' 
a CRR merely creates an intractably large 1ST rate and 
makes it difficult to extract the properties of a CRR from 
those of the composite landscape. Is A=32 large enough? 
The diffusion constants are close to those for |l7| A=108 
and |l(| A = 256, with slightly weaker T-dependence. 
The crystal melts at T m ~ 1.6 an d p| ] fitting D(T) 
yields T c = 0.52, while at 7V=256 |jl6j T m ~ 1.8 and 
(for 'modified' LJ) Q T c = 0.475. The saddle order 
jL2|,[n§ K(T) extrapolates to zero at [[TJ) T c (from D(T) 
by definition), as is so for pi N = 108 and @ for mod- 
ified LJ, N = 256. The plot jlfj of < U ls (T) > is similar 
to that of ref jl3) , reaching the 'bottom' of the landscape 
at the same T w 0.50, and with the gentler drop from the 
high-T plateau expected for small A. Clearly, the essen- 
tial physical features of the crossover - we make no claim 
about deeply supercooled liquids - are present at A=32 
and using a larger system would simply fine-tune vari- 
ous numerical estimates at enormous [ [Tof computational 
cost. 

Natural LJ units will be used throughout. Rather than 
making a single long MD run in the supercooled liquid, 
we average all calculated quantities over an ensemble of 
quenches, starting from different high T — 5.00 configu- 
rations, to avoid (19) broken ergodicity. The hot liquid 
is cooled in one step to a temperature in the 1.20-0.60 
range. The system is equilibrated for 2.5 tlj, data are 
gathered for 62.5 tlj, T is decreased by 0.02 and the pro- 
cess is repeated 10-25 times, generating a single quench 
run; most quenches sampled 16 T. The cooling rate is 
3.08A10~ 4 . Sastry et al. jllj found in an LJ mixture 
that quenches with a cooling rate of 2.70A10~ 4 , close to 
ours, exhibited ~ 75% of the fall in < Ui S (T) > attained 
by quenching almost 100 times slower and reached an 
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apparent bottom somewhat below T c ; this should be ad- 
equate for probing the crossover. At N = 32 the abrupt 
decreases in U signaling solidification, common at N = 
256, do not occur but some quenches develop solid- like 
pair distributions and these are discarded. Results are 
averaged over 23 quenches at the lowest T and 30 at the 
highest. Quench-to-quench fluctuations are much larger 
than any systematic changes over T = 0.02, so we also 
average results at each T with those from the next higher 
and lower T. Even so our data are somewhat noisy but 
the trends are clear for f.10 > T > 0.34, encompass- 
ing the crossover range T a (k, 1.0) > T > T c {= 0.52). 
We believe that our results are well averaged and equi- 
librated down to ~ T c , and are out of equilibrium (but 
'well averaged' over the set of non-equilibrium configura- 
tions allowed by the quench) at the lowest T. 

Conjugate gradient minimizations are performed every 
5 time steps {dt = .00125), or 160 minimizations/TXj. 
Since the range of < u is (T) > is from 9.0 IST/t L j at T 
= 1.10 to 0.23 at T = 0.34, this is sufficient. The equiv- 
alence of conjugate-gradient and steepest-descent for the 
relevant T was verified. The distribution of IST-vector 
(logarithmic) lengths d is bimodal, and [50] to avoid 
counting contributions of two-level systems and anhar- 
monicities irrelevant to diffusion we record a transition 
for d in the large-displacement lobe only. Similar consid- 
erations have [ElJ entered the efforts to relate D to Im—uj 
instantaneous normal modes, where non-diffusive modes 
must be discarded. 

At each T we obtain the three MSD, r2 is (t), R2(t), 
and r2(t) (yielding D), the distribution of waiting times, 
g(t), and the averaged participation ratio < Pr > of the 
IST-vectors. With the i-dependent IS in hand evaluation 
of R2(t) and < Pr > is straightforward. The intrabasin 
MSD r2i s (t) is found from an ordinary MSD algorithm 
with the condition that, when an 1ST is detected, the 
current run is terminated and a new run is begun (new 
origin of coordinates, MSD=0). Binning the times spent 
in an IS, instead of simply calculating the averaged 1ST 
rate, is all that is required for g(t). 

III. LANDSCAPE DYNAMICAL PROPERTIES 

A. Intra- and Inter- Basin Dynamics 

Distributions of IS waiting times g(t) are described 
very well by a KWW (stretched exponential) function 
exp{—(t/r)^) over the entire temperature range. The 
KWW fit is superior to a sum of two exponentials, even 
though the latter has one more adjustable parameter. 
The average (3 = 0.50, and deviations appear to be noise 
with no systematic T-dependence. To investigate trends 
in (3 over a broader T-range we obtained data at normal 
liquid T = 2.00 and also found, to our surprise, j3 = 0.50. 



One hesitates to read very much into the intriguing 
half-integral value of /3 because g(t) is a composite prop- 
erty, dependent upon N crr . An estimate of <?/ oc (i) may 
be obtained as follows. The probability that there is no 
1ST in time t is P no {t) = / t °° dt'g{t'). In a composite, 
no transition means that no CRR has a transition, and 
P no {t) = (/ t °° dt'gi oc (t')) Ncrr . Equating the two expres- 
sions, solving for J" t °° dt' gi oc (t') , and taking d/ dt yields 

,.s ffW h s 

9loc[t) ~ N crr ( f t °° df s (f ))d-V^) • [L) 

The averaged participation ratio of the IST-vectors is 
nearly constant, fluctuating between 19-22 over the T- 
range with an average of 21 for z*=7.0 atoms in a CRR 
and N crr — 4.6 CRR in the simulation box. Using the 
average N crr at all T we calculate gi oc which are also 
well described by a KWW form, with T-independent f3 
= 0.36. Some representative gi oc (t) from Eq ^, simulated 
g(t) and fits with (3 = 0.36 and /3 = 0.50, respectively, are 
presented in Fig [l[ The shift of gi oc (t) to longer times is 




FIG. 1. 'Composite' and local waiting time distribution 
pairs g(t) and gi oc {t) with overlapping KWW fits at T = 0.34, 
0.50 and 1.00, top to bottom. Curves are set to unity at t — 0, 
shifted for display, and decay to zero at the horizontal lines; 
local distribution has slower decay. 

evident. Processes with KWW distributions are not well 
characterized by a single time, but the best compromise 
is the correlation time, n oc = dtgi oc (t)/ gi oc (0). Fig || 
provides further evidence |l(],[l7],[lj| that the low-T dif- 
fusive mechanism sets in at T ~ 0.50, where r/ oc (T) be- 
gins to increase strongly. Also shown is N crr / < loi s >, 
the obvious intensive 'composite corrected' time available 
from the average 1ST rate, in reasonable agreement with 

Tloc{T). 
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FIG. 2. T-dependence of local waiting time from gi oc (+) 
and from N crr / < uii a >(x); smooth curves are spline fits. 



finite size is felt, with r; oc = 0.68 << t p i (thus r2 is (t) 
is noisy at long time). This is illustrated with the inclu- 
sion of gi oc , scaled to intersect r2 is (t) at r/ oc . The 1ST 
are 'bookkeeping' events as border dynamics prevails, D 
is not proportional to < uj is >, there is no effective 
barrier to leaving a basin and hopping is surely not an 
apt description. One might say either that intrabasin dy- 
namics are irrelevant to diffusion, or that no distinction 
exists between inter- and intra- basin dynamics. 

At T = 0.80 < T A , r2 ls (t) looks rather different - the 
difusive interval has grown to the point that it resembles 
a conventional MSD and easily allows estimation of an 
intrabasin Di S . The plateau must be reached eventually, 
but in the current simulation all we can say is t v i > 
6.25. Comparison with r2(t) shows (Figfj]) that D is < D, 
with As=0.0045, D=0.007I4. The two MSD have a brief 
interval of overlap at the beginning of their linear regimes 
but soon diverge, so the mechanism is not governed by 



Let us now view the changes in the mechanism of dif- 
fusion with decreasing T through the paired conventional 
and intra-IS MSD, r2(t) and r2j s (t), combined with our 
knowledge of r; oc (T). The T ~ Ta scenario is found at 
T = f .10, Fig [| a relatively high (although supercooled) 
temperature ~ 2T C where thermal energy is sufficient for 
the system to move freely among closely-spaced IS 




FIG. 3. Conventional (upper) and intrabasin MSD at T 
= f .10. Decaying curve is distribution of local waiting times, 
scaled to cross intra-MSD at t; oc . Plots are closely-spaced 
data points. 

borders. The unconstrained MSD diverges from r2i S (t) 
at t ~ 0.36, before it has even fully established its lin- 
ear diffusive form. For the same reason the trajectories 
which stay in the basin explore it quickly; the short time 
rise in r2j s (t) bends over to form the plateau at about 
T p i ~ 2.0 with only a glimpse of intermediate-t behavior. 
Even so, almost all trajectories leave the basin before its 
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FIG. 4. As Fig 3 for T = 0.80. 

intrabasin diffusion and is probably unchanged from T 
— 1.10. Nonetheless we believe that the emergence of 
intrabasin diffusion is significant and characteristic of the 
upper crossover regime. Basins are not harmonic away 
from the minima and apparently below T4 the higher 
elevations are rough enough to produce diffusion. An 
equivalent statement is that the process is best viewed as 
a random walk among the saddles Jl2|,[i3| . Configurations 
in the upper part of a basin are closer to the saddles than 
to the IS and will be assigned accordingly by the critical 
point mapping. A random walk among saddles connected 
to a basin will appear as intrabasin diffusion. At this T 
ordinary diffusion is also saddle dynamics, so with Di S < 
D it follows that unconstrained saddle dynamics is faster 
than intrabasin saddle dynamics. The need to utilize 
an inefficient random walk to escape the basin, and of 
course the lower thermal energy, has caused the plateau 
time to increase significantly from its T — 1.10 value and 
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i~ioc << Tpi holds even more strongly despite an increase 
in Tioc- 

In our picture of intrabasin diffusion, saddle transi- 
tions occur while the IS does not change. We previously 
discussed 'bookkeeping' 1ST, where the IS changes but 
there is little actual motion. With several IS available 
from a saddle, motion about the basin of attraction of 
a single saddle naturally generates bookkeeping 1ST; the 
IS changes but the saddle does not. Both possibilities, 
which coexist at intermediate T, may be observed with 
a time-dependent critical point mapping, made possible 
by [|l8| a new, very efficient algorithm. Fig || shows IS 



Cooling to T = 0.50 ~ T c produces a remarkable re- 
sult. The period of overlapping linear t-dependence of 
r2 is and r2 is now substantial, and (Fig |^) the true 
diffusion constant may be calculated from trajec- 
tories with no 1ST. The system now stays in a basin 
long enough that on average intrabasin diffusion is fully 
developed before an 1ST occurs, and before the plateau 
is reached (although t; oc is now approaching r p i, Fig||). 
Under these circumstances diffusion is intrabasin domi- 
nated, and the 1ST simply act to keep the process going, 
ie, to avoid the influence of finite basin size. Intrabasin 
and unconstrained saddle dynamics (ordinary diffusion) 
are indistinguishable. We thus have some new 
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FIG. 5. Saddle (upper) and IS energies vs t, T = 0.70. 

and saddle energies vs t at T = 0.70. In the time interval 
0.00 < t < 0.70 the system undergoes seven IS transitions 
while assigned to the same saddle (bookkeeping 1ST) . For 
0.70 < t < 1.15 the IS is constant and there are two 
saddle transitions at t ~ 0.8, followed by rapid exchange 
between a pair of saddles starting just after t = 1.00 
(intrabasin diffusion) . 
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FIG. 6. As Fig 3 for T = 0.50. 

perspective on the changes around T c . There the chain 
of 1ST becomes jn| a Markov process, and an obvious 
explanation is that sufficiently long confinement of the 
system in a basin allows successive 1ST to lose correla- 
tion. While conventionally confinement arises from en- 
ergy barriers, we find that the system diffuses to the bor- 
der and crosses with ease. Confinement is produced by 
diffusion itself; the time required to escape by the ran- 
dom walk is long even absent a barrier at the border. 
A sufficiently slow walking about the saddles connected 
to a basin can produce a Markov process among the IS 
themselves, and IS dynamics is also a good description 
of the intrabasin dominated regime. Upon close exami- 
nation, the crossover is not simply from saddle dynamics 
at high T to IS dynamics at low T. There is a third, 
intermediate-T region where both descriptions are cor- 
rect, defined by the condition D — Di s . 

Arrhenius plots of D is and D are displayed in Fig 0. 
Estimates of Di 8 are possible for 1.10 > T > 0.38. Lin- 
ear behavior of r2i S (t) is not cleanly visible below T = 
0.38 and is marginal at T = 1.10; we expect it disappears 
at slightly higher T, but collecting data where trajecto- 
ries leave a basin so quickly is difficult. The intrabasin 



G 



dominated D = Dj S range, marked by arrows, is approx- 
imately T c = 0.52 > T > 0.46. Starting from high T, the 
upper limit of intrabasin dominance is reached when two 
conditions hold: T is low enough for r2,* s (i) to have a 
diffusive ^interval, and t; oc adequately exceeds the time 
at which that interval begins. For the lower limit, ei- 
ther of two conditions may be true: r/ oc substantially 
exceeds t p i, so a barrier to leaving the basin exists and 
activation becomes important, or the linear portion of 
r2i s (t) shrinks so that it no longer represents the intra- 
basin motion. In fact, both appear to set in gradually 
and simultaneously. 




1 1.4 1.8 2.2 2.6 
1/T 

FIG. 7. Natural logarithms of intra-basin diffusion con- 
stant Di S (dashed line spline fit and points) and true D (solid 
line) vs 1/T. Intrabasin dominated D = D,: s range delimited 
by arrows, left arrow coincides with T c . 

At our lowest temperature, T = 0.34 (Fig ||), r2 is (t) 




rises quickly to the plateau just as at the high T = 1.10; 
a linear region of this curve cannot be identified and t p i 
is now back within the 6.25 tlj window. 

Quantitative estimates of t p i(T), along with r; oc (T), 
are shown in Fig O; at intermediate T all we can say is 




FIG. 9. Local waiting time r; oc and, in two segments (line 
plus ★), rough estimate of plateau time t v i vs T; low-T data 
shown with quadratic fit. Arrow indicates T c . 

T p i > 6.25. While t p i is short at high T because of the 
abundant thermal energy, it is short at low T because 
the system occupies the lower regions of the basin only. 
Thus a) the plateau is reached after a shorter displace- 
ment, and b) the intra-basin motion is mostly fast har- 
monic oscillations, not diffusion. The t p i and r/ oc curves 
cross at T w 0.51, and we propose r p i — ti oc as a new 
way to define a crossover temperature. Below T = 0.51 
the system explores the accessible portion of the basin 
more quickly but exits more slowly, indicating barriers 
at the border. Gradually saddle dynamics loses out to 
IS dynamics as the physical description, D ^ D is below 
T = 0.46, and barrier crossings assume the rate-limiting 
role for diffusion. 

In addition to the standard analysis of D(T) and the 
results |l3|,|l(],|l7| of very recent work, we now have two 
new indicators pointing to a change in the mechanism of 
diffusion at T ~ 0.50. The local waiting time r/ oc crosses 
t p i at T = 0.51, and intrabasin dominance sets in at T 
= 0.52 with D = D ls for 0.52 > T > 0.46. However 
the crossover develops over a wide temperature range, 
and involves much more than the presence or absence of 
hopping, as explained in the last few paragraphs. We 
hope it is clear that the IS-mapping is a powerful tool for 
studying the dynamical crossover. 



FIG. 8. As Fig 3 for T = 0.34. 
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B. Landscape Vector Correlations, Saddle Dynamics 
and Borderism 



IST-vectors are anticorrelated |L0|] at high T, and cor- 
relation decreases through the crossover range until the 
IST-Markov approximation becomes accurate below T c . 
A substantial negative value of C{[3) is a sign of saddle 
or border dynamics, while C{0) ~ indicates IS or hop- 
ping dynamics. IS borders are closely spaced near a mul- 
tidimensional saddle. Anticorrelation is required for the 
'bookkeeping' 1ST arising from a small motion through 
such a region to not, wrongly, predict a large motion. It 
implies that saddle dynamics is a more physical descrip- 
tion than IS dynamics. A simple process appears com- 
plicated when formulated with inappropriate elementary 
steps, and loss of correlation signals the onset of the low- 
T mechanism, which is truly IS dynamics. 

Analysis of the MSD within the basin of attraction of 
an IS, r2 is (t), proved quite fruitful. We will now see 
that such is also the case for the MSD of the IS itself, 
R2(t), leading to another informative correlation, that 
of the return vector and the IS configuration. With the 
definitions from the Introduction, 



r2(t) = R2(t)+ < (Au(t) f > /6N + 
< AR(t) • Au(t) > /6N, 



(2) 



It follows that r2(t) > R2{t) unless the IS displacement, 
AR(t), and the return vector displacement, Au(i), are 
anticorrelated. Fig |o| shows that at T = 1.00, R2(t) lies 
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FIG. 10. Ordinary MSD (solid) and IS-MSD (dashed) vs t; 
top to bottom, T = 1.00, 0.90, 0.80, 0.70. 

above r2(t), they nearly sumerimpose at T = 0.80, and 
the 'obvious' r2(t) > R2{t) holds at T < 0.80. Thus the 
pattern of anticorrelation of IST-vectors <5R a • <5R Q+j a at 
high T, decreasing with decreasing T, is repeated with 
the AR(i) • Au(i) correlation. 



For a more quantitative analysis we have calculated 

< (Au(i)) 2 >. One anticipates a rapid decay of < 
u(t) • u(0) > leading to < (Au(t)) 2 >^ 2 < q 2 > and 
this expectation is correct. The correlation Cjt u (i) =< 
AR(t) • Au(f) > /6N, determined from r2(i), R2(t), and 

< (Au(t)) 2 > via Eq |[ decays quickly (Fig [ll]) to an 
asymptotic negative value which decreases in amplitude 
with decreasing T. We suggest that Cr u is a direct mea- 
sure of the degree of borderism - the extent to which the 
system 'lives' on the IS-borders [E2|. 
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FIG. 11. Correlation Cn u (t), amplitude increasing with in- 
creasing T; top to bottom, T = 0.60, 0.70, 0.80, 0.90, 1.00, 
1.20. 

Expressing AR(t) as the sum of the IST-vectors the 
most important contributions are evident, 



CruW =< SR n ■ u(t) > - < 5Ki u(0) > 
+ (smaller terms). 



(3) 



The first term on the RHS may be understood with 
Fig [l2| a one-dimensional illustration of the situation at 
time t, just after the n'th 1ST. The system has moved 
from the left- to the right-hand basin and it follows that 
<5R„ and u(t) point in opposite directions, with a neg- 
ative contribution to Cr u . For the second term on the 
RHS, set t = in Fig The first 1ST will occur from 
right to left in the future so change 5R„ to <5Ri and 
reverse its direction. The vectors now point in the same 
direction so with the minus sign in Eq|^ this contribution 
to Cr u is also negative. In the 3N dimensional config- 
uration space (or 3N crr dimensional local region space) 
the vectors will be distributed about the d = 1 relative 
orientations, reducing the dot products. 
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FIG. 12. Return vector u(t) and IST-vector <5R„ at time 
t just after n'th (most recent) 1ST. 

The geometrical arguments apply when u(0) is evalu- 
ated just before, and u(i) just after, an 1ST. For border 
dynamics these conditions 'always' hold and Cr u < 0. 
On the other hand if the system oscillates about the IS 
for long periods between 1ST, the low-T scenario, the re- 
turn vector has no special relation to any 1ST vector at 
any time and Cr u ~ 0. The correlation of AR(t) and 
Au(i) is a maximum for pure border dynamics and van- 
ishes in a hopping model. The interesting part of Cr. u is 
its asymptotic value and we define 

B(T) = -C Ru (oo,T) (4) 

as the border index. The T-dependence of the border 
index is shown in Fig It extrapolates to zero via 
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FIG. 13. Border index B(T) (points) and power law fit, 
tending to zero at T — 0.53. 

power-law fit for 1.2 > T > 0.60 at T = 0.53, yielding 



another physically appealing route to a crossover tem- 
perature. The fit range differs from that used for T c but 
B(T) is close enough (power = 0.74) to linear in T that 
we expect only a small effect. 

IV. DISCUSSION 

In the foregoing we have examined the T-dependent 
physical mechanism of diffusion in unit-density LJ. For 
some perspective, it is useful to ask why this is a non- 
trivial problem. LJ is classified as weak/nonfragile, 
indicating weak T-dependence at high T and a constant 
activation energy at low T . The simulated diffusion con- 
stant is indeed well represented by 

D(T) = 0.39Tea;p(-1.16/T) (5) 

for 1.10 > T > 0.34; E A = 1.16. The factor of T is 
appropriate because D ~ T at high T and constant den- 
sity. Linear or other multiplicative powers are not very 
important when fitting over several decades but we have 
only 20X and they significantly influence the goodness of 
fit and the value obtained for the activation energy. De- 
spite the simplicity of Eq ||, fitting |lj| the same data to a 
power law for 2.0 >T> 0.60 yields a reasonable [|l3| T c . 
One might well wonder if anything is going on beyond 
a gradual change from T » Ea to T « Ea- That 
is a bare-bones prescription for a 'transition to hopping', 
consistent with the 1ST approaching a Markov chain JjC[ ] 
and the saddle order extrapolating to zero [ fj^]l8| . 

However it then seems odd that Ea — 1.16 does not 
appear as a characteristic dynamical temperature instead 
of various T ~ 0.50, and the sharp drop in < Ui S (T) > 
around T c suggests that the crossover is more interest- 
ing. To investigate this crucial point we have begun |lq] 
to catalog the saddles and barrier heights in unit-density 
LJ. Starting from the thermal configuration the IS is de- 
termined and the connected first-order saddles are found 
fl8f by eigenmode following. Steepest descent from the 
saddle yields barrier heights AU, which may be collected 
as a function of IS energy |23| ] or of the original tempera- 
ture. Preliminary results Jl8| are that the barrier heights 
increase through the slope region |||ll| of < Ui S (T) > 
from < AU (T = 1.20) > = 2.8 to < AU(T = 0.40) > 
= 3.7; T is less than the averaged barrier height for the 
entire range of this study. The same result has recently 
been obtained |23| in a soft-sphere mixture. Comparisons 
of barrier heights to T must be made carefully. Most 
discussions of barrier crossing assume that AU is much 
greater then T, but that does not apply here and one 
might argue that the relevant height is (AU — T). Even 
with such a modification, average barrier heights remain 
larger than Ea, indicating - along with the new results 
found above - that the mechanism of diffusion in LJ is 
highly non-trivial. 
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There are two explanations for the clear presence of 
a high-T mechanism with free motion among the basins 
at T > T c despite an apparent requirement for activated 
hopping. One is the broad distribution of barrier heights, 
reflected in the low (3 = 0.36 for the KWW distribution 
of local waiting times. Even if T falls below < AU(T) >, 
there will be no transition to low-T dynamics if enough 
low barriers remain. The average alone does not ade- 
quately describe the physics. Furthermore there are sev- 
eral ways to perform the average pl| and a physically 
motivated choice must be made. For example a calcula- 
tion of < AU(T) > including saddles of all orders would 
be nonsense, since high barriers which are never visited in 
thermal motion would dominate. Averaging over first or- 
der saddles eliminates many irrelevant high barriers, but 
not necessarily all. Other topological features such as the 
connectivity of the IS and saddles should be important 
as well. The landscape determines T c and more generally 
D(T) through the interplay of several physically signif- 
icant effects. Note that one important landscape quan- 
tity < Ui S (T) >, is a 'thermodynamic' property of the 
basins. One can imagine different liquids with the same 

< Ui S (T) > and different T-dependent distributions of 
barrier heights, giving rise to quite different activation 
energies and other aspects of dynamics. 

Another is found in Cavagna's |l2|,^3| argument that 
barrriers are irrelevant while the saddle mechanism dom- 
inates, so the action of a high-T mechanism while T is 
less than < AU > need present no conundrum. Then, 
when the saddle mechanism shuts down and activation 
becomes required at T ~ 0.50, the barrier heights have 
already completed most of their growth. The system 
will exhibit approximate Arrhenius behavior upon fur- 
ther cooling. Cavagna gives a 'fragile' scenario where 
the barriers are already >> T at T c , for a very strong 
increase in relaxation time. However although we have 

< AU{T C ) > /T c ~ 6, E A is only - 2T C , presumable due 
to the importance of low-barrier paths. Relatively weak 
Arrhenius behavior at T < T c results, consistent with LJ 
being [^) a weak/nonfragile liquid. 

Even if Ea — 1-16 is not meaningful at higher T, Eq|| 
can still fit the data because the Arrhenius factor is ap- 
proaching its ^^-independent high-T limit. Saddle dy- 
namics, in which D(T) is pif| proportional to the number 
of Im — u> instantaneous normal modes or [fl7[ the saddle 
order < K(T) >, apparently is consistent with D ~ T. 
Thus Eq U can represent the full T-range, but the empir- 
ical activation energy is a compromise resulting from the 
fitting process with no simple relation to a barrier height. 

In either case the crossover is seen to be a subtle, com- 
plex process. Application of the IS-mapping has revealed 
far more then could be deduced from D(T) alone, includ- 
ing: 1. intrabasin and intrabasin dominated diffusion, 
associated with saddle dynamics 2. the distribution of 
local waiting times, well behaved as N — > oo, calculated 
from the simulated distribution with the idea that the 



configuration space is a composite 3. a new, physically 
appealing criterion for dynamical crossover, r; oc = t p i 4. 
a quantitative indicator B(T) of the dominance of sad- 
dle or border dynamics (borderism), which also yields a 
crossover temperature. 

The value of T c = 0.52, from a power-law fit to D(T), 
is now just one of many pieces of evidence that a dy- 
namical crossover occurs at T ~ 0.50. From recent prior 
work |l0|]i~7| ]: 1. IS energy < Ui S (T) > is undergoing 
its steepest fall. The true bottom is uncertain due to 
cooling-rate dependence and point 1. only identifies the 
general vicinity of the crossover. 2. saddle order K(T) 
also extrapolates to zero at T — 0.52 3. IST-Markov 
approximation becomes accurate at T ~ T c . In this pa- 
per we have added: 4. D = Di S , diffusion is intrabasin 
dominated, for 0.52 > T > 0.46. 5. local waiting time 
reaches time to explore basin, r/ oc =t p i, at T — 0.51. 6. 
border index B{T) extrapolates to zero, T — 0.53. We 
suggest that points 3-5 are particularly direct evidence 
that the mechanism of diffusion is changing. They are 
based on ideas about the mechanism itself, while a fit to 
D(T) produces a characteristic temperature based upon 
the 'symptoms'. Although 1-6 focus attention on a nar- 
row range 0.53 > T > 0.46 the crossover is gradual and 
might be thought to take place over 1.10 > T > 0.38, the 
interval where intrabasin diffusion can be detected, en- 
compassing the slope of < Ui S (T) >. This agrees roughly 
with having the crossover begin at T A - 

We believe it significant that intrabasin diffusion can 
produce confinement for times long enough that the 1ST 
become a Markov chain, without invoking barriers. Thus 
diffusive confinement and confinement by barriers be- 
come two distinct features of the mechanism. Since in- 
trabasin diffusion can be explained as saddle dynamics, 
our work fits with the growing recognition that the sad- 
dles of the landscape need to be considered on an even 
footing with the IS. There is no contradiction in having 
the intrabasin diffusion and IST-Markov regimes over- 
lap. In the latter case D oc< uJi S >, and in the former, 
since escape from a basin requires diffusion to the border, 

< LO is >cx D is . 

Our simulation model exhibits most of the physical 
features currently being discussed for moderately super- 
cooled dynamics, as well as some found here for the first 
time. There is much to be learned from LJ, N = 32, 
although of course we intend to apply our methods to 
other liquids. It is difficult to achieve equilibrium below 
T c . Note that, even though we report data down to T 
= 0.34, only one significant T appearing in the above es- 
timates is below 0.50, the end of intrabasin dominance 
at T = 0.46. Our interest is the crossover and the low- 
est, probably non-equilibrated, T are not required for the 
principal conslusions. 
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